Total energy density as an interpretative tool 
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We present an unambiguous formulation for the total energy density within density-functional 
theory. We propose that it be used as a tool for the interpretation of computed energy and electronic 
structure changes during structural transformations and chemical reactions, augmenting the present 
use of electron density changes and changes in the Kohn-Sham local density of states and Kohn-Sham 
energy density. 
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Changes in the electron density of a system undergo- 
ing a structural transformation or a chemical reaction 
have long been used to understand the driving forces 
underlying the transformation or the reaction. Indeed, 
the electron density itself has been regarded as one of 
the most significant descriptors of a system of electrons 
and nuclei since the advent of quantum mechanics. This 
view of the importance of the electron density was greatly 
strengthened by the emergence of density-functional the- 
ory. The Hohenberg-Kohn theorem [jjj showed that the 
total energy of a system could be treated as a func- 
tional of the electron density. The Kohn-Sham theory 
H makes possible the construction of the true ground- 
state electron density of an iV-electron system from the 
N one-electron orbitals of lowest energy of a particu- 
lar independent-electron system. This rigorous decom- 
position of the electron density into contributions from 
one-electron orbitals, the Kohn-Sham orbitals, greatly in- 
creases its power as an interpretive tool. 

At present, methods based on density- functional the- 
ory are dominant in first-principles calculations for con- 
densed matter and for much of quantum chemistry. The 
electron density and its Kohn-Sham orbital decomposi- 
tion, fundamental outputs of the computations, are typ- 
ically used to enrich understanding of the resulting ener- 
getics. For example, studies of the isosurfaces of electron 
density changes along a reaction pathway for the disso- 
ciation of H2 on the Pd (100) surface || have revealed 
the formation of bridge bonds between the s-p tails of 
the metallic surface electron density and the a g and a* 
molecular orbitals. These bridge bonds mediate the hy- 
bridization of the molecular orbitals with the metallic 
d-orbitals before they actually overlap. It is the bridge 
bonds which evolve into the bonds between the dissoci- 



ated hydrogen atoms and the metal. 

The Kohn-Sham theory enables one to go still deeper 
into the interpretation of the results of the electronic- 
structure computations. For extended systems, one can 
construct the Kohn-Sham density of states, 
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where ej is one of the N lowest Kohn-Sham eigenvalues, 
and use it to tease out the energetics associated with the 
electron-density changes. For example, in studies of the 
interaction of chemisorbed O and H to form OH and H 2 
on the (111) surfaces of Rh and Pt the resonances 
and bound states associated with the bonding and anti- 
bonding orbitals of the atoms, radicals, or molecules to 
the surface are shown to be clearly visible in g(e). Even 
more revealing is the local density of states, 
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where 4>i(v) is one of the occupied Kohn-Sham orbitals, 
which enables one to associate a particular resonance or 
bound state in the Kohn-Sham spectrum, Eq.(Q), with 
a particular atom or molecule. It allows one to make 
a resolution with respect to the Kohn-Sham energy of 
the changes of electron density accompanying structural 
transitions or chemical reactions. 

Nevertheless, what drives a transformation or a reac- 
tion is the dependence of the total energy on the nuclear 
coordinates. Thus, what one needs as an additional inter- 
pretive tool is an unambiguous, computationally feasible 
total-energy density e(r), such that 
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E = / d 3 re(r), 



(3) 



where E is the total energy. With e(r) one could explore 
the relation between local changes in the electron density 
and the spatial dispersion of the total energy changes 
through the corresponding changes in e(r). 

To be more specific, in Ref. |J, evidence was found 
for the simultaneous presence of all the specific mecha- 
nisms commonly cited as playing important roles in the 
dissociation of H2 on a transition metal surface. These in- 
cluded the formation of an orthogonality hole in the s — p 
tail of the metallic electron density, with concommitant 
flow of screening charge into the d-states, as emphasized 
by Harris and Andersson j^] ; the hybridization of the te- 
states with the occupied a g bonding orbital of the H2, 
as emphasized by Hammer and Scheffier Q, and the hy- 
bridization of the unoccupied a* antibonding orbital of 
the H2 with associated backbonding, as emphasized by 
Hammer and N0rskov |t]|| . Each of these processes intro- 
duced distinct, easily recognized characteristic changes 
in the electron density, which occured in distinct regions 
of space. One could therefore utilize e(r) to order the 
relative importance of each of these mechanisms as well 
as of bridge-bond formation at different locations along 
the dissociation pathway, garnering thereby a detailed, 
intimate quantitative understanding of the breaking of 
the intramolecular bond and the formation of metal-atom 
bonds in the course of the dissociation. 

There are several difficulties to be surmounted in con- 
structing an unambiguous, computationally- feasible ex- 
pression for e(r), as will be discussed below. To avoid the 
difficulties, only the Kohn-Sham energy density e s (r), 
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was used as an interpretive tool in reference (J], though 
to good effect. In the present paper, we present an ex- 
plicit formulation for e(r) which is both unambiguous 
and computationally feasible, as required. 

The total energy can be represented as a density func- 
tional in the form 



E = T + V et 



Vn 



(5) 



where T, V ee , and V en are the functionals for the elec- 
tron kinetic energy, the electron-electron interaction, and 
the electron-nuclear interaction. V nn is the Coulomb in- 
teraction between the nuclei. In Kohn-Sham density- 
functional theory ]9|Jl0|| , the expression Eq. (||) is rewrit- 
ten as 



E = T s + U + E xc + V en + V n - 



(6) 



Here T s is the kinetic energy of the independent Kohn- 
Sham particles, U is the Hartree electrostatic energy, 



U= l - [d*r fd\' n ^ n{r [\ 
2 J J r-r' ' 



(J) 



where n(r) is the electron density, and E xc is the 
exchange-correlation energy, 



E xc = T - T s + V cc - U. 



(8) 



On the other hand, the total Kohn-Sham energy of the 
system is 



N 

E S =J2^=T 3 + J d 3 rn(r)v s (r), 



(9) 



In Eq. (||) , v s (r) is the Kohn-Sham potential, 

u s (r) = v en (r) + u H (r) + v xc (r), (10) 

where v en (r) is the Coulomb potential produced by the 
nuclei at r, which enters V m as well 



V en = / d 3 rn(r) v en (r), 



(11) 



v H (r) is the Hartree potential, the mean electrostatic po- 
tential produced by the electrons, 
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which enters the Hartree energy functional, 



(12) 



U = J Ae H (r), e H (r) = in(r)w H (r). (13) 

Finally, ^ X c( r ) is the functional derivative of E xc : 

SE XC 



^xc(r) 



Sn(r) 



(14) 



Thus Eq.@ can be transformed into 

E s =T S + 2U + V en + [ d 3 rn(r)v xc {r). (15) 



Inserting Eq.([l5|) into Eq.(||) allows us to eliminate T s , 
resulting in 

E = E s - U + V nn + E xc - J d 3 rn(r) v xc (r). (16) 

Each of the above terms can unambiguously be written 
as an integral over an energy density so that 



e(r) = e s (r) - e H (r) - ^p„(r) v en (r) 
+ e xc (r) - n(r) v xc (r). 



(17) 



In this expression, e s (r) is given by Eq.(^) and can read- 
ily be constructed from the output of the standard codes 
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which solve the Kohn-Sham equations. The electron den- 
sity n(r) is a fundamental output of such computations. 
The potentials v H (r) and v en {r) must be constructed dur- 
ing the computations, as must be i? xc (r). The number 
density of the nuclei weighted by their respective atomic 
numbers is 



(18) 



where R Q is the position and Z a the charge of nucleus a. 
We note that the third term in Eq. (|l7|) , though singular 
at the nuclei, vanishes between the nuclei. 

The exchange-correlation energy density e X c( r ) is de- 
fined in terms of v xc [n](r) via the procedure of Burke, 
Cruz, and Lam (BCL) pd|], a n extension of the original 



ideas of Engel and Vosko 

V 2 e xc (r) = 3V(n(r) \7v x 
where 
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and n 7 (r) = 7 3 n(7r). Their definition has several advan- 
tages in this context. Most importantly, the BCL proce- 
dure defines an exchange-correlation energy density un- 
ambiguously in terms of any given exchange-correlation 
energy functional (including the exact one), since it uses 
the potential. Thus exchange-correlation energy densi- 
ties calculated this way are approximations to an ex- 
act quantity. Within the local density approximation 
(LDA), their procedure reproduces simply e x " lf (n(r)), 
the exchange-correlation energy density of a uniform elec- 
tron gas of density n(r). This is just the conventional 
energy density within LDA, and so is already calculated 
in standard LDA Kohn-Sham calculations. Within any 
generalized gradient approximation (GGA) such as PBE 
fyl or BLYP their procedure defines a different 

energy density from the conventional one. This energy 
density includes dependencies on the Laplacian and other 
higher derivatives of n(r), and so is far more sensitive to 
details in n(r) than the conventional GGA forms. The 
scaling defined in Eq. (|2(i| ) is easily performed on any 
approximate density functional. Finally, the BCL proce- 
dure defines an exact unambiguous exchange energy den- 
sity, so that it can be applied even to hybrids of GGA 
with exact exchange ]Tl| ]. 

To illustrate this energy density, we plot the various 
contributions to the total energy density for the Na atom. 
Our calculations are all for exact exchange only, using the 
atomic optimized effective potential (OEP) code of En- 
gel |16|. From Figs. 1 and 2, we see that the total 
energy density is usually dominated by two terms: the 
Kohn-Sham eigenvalue contribution and the electrostatic 
energy correction. The two exchange terms of Eq. (|17|) 
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FIG. 1. Radial contributions to the total electronic energy 
density for the Na atom (atomic units). 
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FIG. 2. Radial exchange contributions to the total elec- 
tronic energy density for the Na atom (atomic units); the 
dashed curve is the sum of the other two. 
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FIG. 3. Same as Fig. 1, but for the difference between 
the neutral atom, Na, and the ion, Na + ; the dashed curve is 
Ae(r) - Ae H (r). 



are small, and largely cancel, as shown in Fig. 2. This 
cancellation is already apparent in LDA exchange, where 
the sum of these terms is only 1/3 of ex mf - This agrees 
with arguments [l7-19| that changes in the energy can 
be largely understood in terms of changes in KS eigenval- 
ues and the electrostatic potentials. The present analysis 
opens the possibility of extending these arguments to the 
exact functional. 

To illustrate the importance of the electrostatic energy 
correction, consider Fig. 3, which shows energy density 
differences between neutral Na and its ion. The Kohn- 
Sham contribution is dominated by the core (r < 2, 
roughly). This is because the 3s electron in Na induces 
an almost constant shift in the core Hartree potential, so 
that all eigenvalues are about 0.29 higher in the neutral 
relative to the ion. When the Hartree correction is made, 
and the total energy density difference plotted, we find 
that there are two almost equal and opposite contribu- 
tions to the energy difference. To see that these are en- 
tirely electrostatic effects, we further subtract e H (r) from 
the total (dashed line), showing that, in the absence of 
Hartree contributions, the total energy density lies al- 
most entirely in the valence region. This simple analysis 
demonstrates both the importance of electrostatic contri- 
butions to the energy density, and how the various terms 
are needed in Eq. (|]j]) to produce the final physical pic- 
ture. (For example, overall shifts in KS eigenvalues do 
not contribute to e(r) , but do show up in e s (r) .) The elec- 
trostatic energy correction is omitted in the generalized 
perturbation method widely used in the theory of 
alloys (2|] . The present analysis allows us to explore the 
role of exhange-correlation energy differences in chemical 
systems. 

Note that v H (r) and v en (r) each diverge separately in 
the thermodynamic limit for extended systems. The di- 
vergences cancel each other in v s (r), Eq.(O), because 



only their sum enters. The divergences similarly cancel 
in E, Eqs.(|), @, and ^ and in E s , Eq.(g). They do 
not cancel in e(r) as written in Eq. (|l7|) . This difficulty 
should be obviated by the methods of eliminating the 
singular parts of v u (r) and v en (r) which are standard in 
the theory of the electron structure of extended systems. 
That is, v H (r) and v en (r) should be interpreted as the 
nonsingular part of the Hartree and nuclear potential, 
respectively, in Eq.([l7]) for extended systems. 

In conclusion, we propose (1.) that the total energy 
density between two nuclei be constructed from the out- 
put of Kohn-Sham computations augmented by the pro- 
cedure of ref. |]llf and used as an interpretive tool in an- 
alyzing energy and electronic structure changes during 
transformations and reactions when the focus is on the 
study of bond formation in the spaces between atoms. 
We propose in addition (2.) that the full expression 
for e(r), Eq.(|l7|), be coarse-grained by integration over a 
small sphere containing a particular nucleus to examine 
the role of the corresponding atom in chemical bonding 
or structural transformation. Such spheres are already 
present in, e.g., linear-augmented-plane-wave, muffin-tin, 
and KKR codes and can be readily introduced into other 
procedures. 
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